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Size matters: the non-universal density profile of subhaloes in SPH 
simulations and implications for the Milky Way's dSphs 
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ABSTRACT 

We use dark matter only and full hydrodynamical Constrained Local UniversE Simulations 
(CLUES) of the formation of the Local Group to study the density profile of subhaloes of 
the simulated Milky Way and Andromeda galaxies. We show that the Einasto model provides 
the best description of the subhaloes' density profile, as opposed to the more commonly used 
NFW profile or any generalisation of it. We further find that the Einasto shape parameter 
tie is strongly correlated with the total subhalo mass, pointing towards the notion of a non- 
universality of the subhaloes' density profile. We observe that the effect of mass loss due to 
tidal stripping, in both the dark matter only and the hydrodynamical run, is the reduction of 
the shape parameter between the infall and the present time. Assuming now that the dSphs 
of our Galaxy follow the Einasto profile and using the maximum and minimum values of 
from our hydrodynamical simulation as a gauge, we can improve the observational constraints 
on the R m ax-V max pairs obtained for the brightest satellite galaxies of the Milky Way. When 
considering only the subhaloes with —13.2 <; My <J —8.8, i.e. the range of luminosity of the 
classical dwarfs, we find that all our simulated objects are consistent with the observed dSphs 
if their haloes follow the Einasto model with 1.6 <^ He <^ 5.3. The numerically motivated 
Einasto profile for the observed dSphs will alleviate the recently presented "massive failures" 
problem. 

Key words: methods: numerical - iV-body simulations - galaxies: formation - haloes - Local 
Group 



1 INTRODUCTION 

While the predictions of the current A Cold Dark Matter (ACDM) 
model have been widely confirmed at cosmological scales, there 
are still a number of discrepancies between theory and observations 
at galactic and subgalactic scales: one example is th e well-known 
"missi ng s atellite prob l em", first pointed out by iKlypin et al.l 
d 19991) and lMoore etall jl999h . The high number of substructures 
resolved within the virial radius of galaxy-type objects in high reso- 
lution cosmological simulation mismatches the number of observed 
satellite galaxies of our Milky Way (MW) and nearby galaxies. To 
alleviate the problem one must invoke some mechanisms, such as 
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early reionization of the intergalactic medium and supernovae fe ed- 
back jBullock et alj200d : ISomervillel2002l ; lBenson et alj2002h . to 
suppress galaxy formation below a certain scale. 

However, there is an inconsistency not only with the num- 
ber, but also about the kinematics of the observed MW's dwarf 
spheroidals (dSphs) when compared to the velocity profiles of 
the most massive subhaloe s found in dark matter simulations 
jBovlan-Kolchin et alj[2oTlh . Assuming that these subhaloes fol- 



low a Navarro et al. ( 19961 NFW herafter) profile, they have been 



found to be too dense to host th e MW's bright sa t ellites. This is di- 
rectly related to the findings of iBovill & RicottH J201 lal fbh. whose 
simulations showed an overabundance of bright dwarf satellites 
(L v > W 4 L sun ) with respect to the MW's dSphs. 

A number of studies tried to reconcile simulations with obser- 
vations. 

The po ssibility that the MW is a s tatistical outlier has been 
ruled out by IStrigari & Wechslerl d2012h . who used data from the 
Sloan Digital Sky Survey to show that, down to the scale of Sagit- 
tarius dwarf, our Galaxy is not anomalous in its number of classical 
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satellites. Further, the analysis of iBovlan-Kolchin et alj l l2012h . in- 
dependent from the choice of actual density profile of the simulated 
subhaloes, demonstrates that supernova feedback is unlikely to be 
an explanation for the low inferred densities of dSphs, due to their 
small stellar masses. Different hypothesis for the nature of dark 
matter can naturally lead to the format ion of less concent rated sub- 
haloes in a warm dark matter scenario dLovell et al.l2012l) or in sim- 
ulations of self-interacting dark matter models (IVogelsberger et al.l 
providing an interesting alternative to the ACDM model. 
Moreover, the discrepancy between observed and simulated satel- 
lite galaxies may reflect the fact that the MW is less massive than is 
commonly thought: a total mass between 8 ■ 10 11 <; M/M© <; 
10 12 h as been argued in IVera-Ciro et al.l d2013l) and IWang et ail 
j2012l) . However, lowering the mass of the MW still do not explain 
why its dSphs (as well as many isolated dwarf galaxies), seem to 
live in haloes whose mass is smal ler than the current expectation 
from abundance matching models dFerrero et al.ll2012h . 

The inclusion of baryons in simulations has also been ex- 
plored, and it has been found either to h ave negligible eff ects 
on the dark matter density of subhaloes ( IParrv et al.l 1201 2|) or 
to ha ve a twofold effect on their density profile ( Di Cintio et al .1 
I20T 1). In fact, recognizing that at galactic scales baryonic pro- 
ce sses are expected to p lay a crucial role, it has been investigated 
in lDi Cintio et al. l|) the effect of the inclusion of baryons in 
SPH simulation within the CLUES projectQ These simulations are 
designed and constrained, respectively, to reproduce as closely as 
possible the actual observed Local Group with its two prime galax- 
ies MW and Andromeda (hereafter also referred to as M31) and 
hence serve as an ideal testbed for investigating the dynamics and 
kinematics of the satellite populations of the real MW and M3 1 . In 
this previous study it has been found that, while in some cases the 
baryons are able to lower the central density of subhaloes, through 
mechanisms such as gas outflows driven by star form ation and 
supernovae jNavarro et aill 19961. iGovernato etaTll2012l) . there are 
still substructures whose densi ty is increased, as expect ed from the 
adiabatic contraction model of lBlumenthal et al (1986). 

The underlying assumption in many previous works is that 
the satellite galaxies of the MW are embedded in subhaloes whose 
mass profile is described by the NFW model: it is still a matter of 
debate, however, if this profile is the best choice in modeling the 
dSphs' density. 

On one hand, IWalker & Penarrubial d201ll) constructed a 
method for measuring the slope of the mass profiles within dSphs 
directly from stellar spectroscopic data, independently from any 
dark matter halo model and velocity anisotropy of the stellar trac- 
ers, and showed evidence for the profile of the Fornax and Sculptor 
dSphs to be consistent with cores of constant density within the 
central few-hundred parsecs of each galaxy, thus r uling out a cuspy 

[rofile such as the NFW one. On the other hand, IWolf &Bullockl 
20121) used a Jeans analysis to show that, even in the limiting case 
of an isotropic velocity dispersion, not all of the dwarfs prefer to 
live in halos that have constant density cores . It mu st be noticed, 
however, that while the lWalker & Penarrubial d201 lb method is in- 
sensitive to orbi t al ani sotropy and underlying halo potential, the 
IWolf & Bullockl J2012h results are dependent from these yet un- 
known quantities. 

Regarding simulations. |Pi Cintio et al. I J20T l]) pointed out that 
the NFW profile may not be appropriate to describe the sub- 
haloes' density. While computing the subhaloes' circular velocities 
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IBovlan-Kolchin et all f20ll) used an Einasto profile to model the 
density distribution of subhaloes at small radii, in order to correct 
for the effects of the force softening, and the raw particle data at 

higher radii. 

The three-parameter Einasto profile ((Einasto 1965|), indeed, 
has been found to more accurat ely describe the halo density in 
dark matter only simulat i ons (e.g. [ Navarro et al. 2004; Merr itt et al.l 
20061: IPrada etai1l2006l: iGao et al. 2008; Havashi & White! 120081: 



Navarro et al.ll2010l; iLudlow et al.ll2011l ; iReed et al.ll2011l) . reduc- 



ing the residuals of the fits by 20% with respect to the correspond- 
ing NFW models. 

In this Paper we mainly focus on hydrodynamical simulations 
and, after a brief explanation, in Section [2] of the CLUES project 
simulations, we study the mass profile of substructures within the 
two main haloes of the simulated Local Group, formally calling 
them Milky Way and M31. In Section[3]we focus on the quality of 
several analytical models in describing the density profile of galac- 
tic subhaloes showing that, also in hydrodynamical simulations, the 
Einasto profile provides the best description. In Section|4]we then 
show that the corresponding profile shape parameter ?ie scales with 
the virial mass of the subhalo. We finally discuss the implications 
for the mismatch between the kinematics of the observed MW's 
dSphs and the simulated substructures in Section [5] before con- 
cluding in Section[6] 



2 THE SIMULATIONS 

Our simulations form part of the aforementioend CLUES project 
and are based upon a WMAP3 cosmology. These constrained 
simulations of the Local Universe have already been presented 
and extensively used for other investigations an d we refer the 
reade r to those articles for more d etails (e.g. iGottlober et all 



20101: iLibeskind et al.l |2010L l201ll ; iKnebe et all |201C| . 



2011; 



Di Cintio etalj|20 lib . We therefore only repeat here the most ba- 



sic informations. Th e simulations assume a WMAP3 cosmology 
dSpergel et aill2007l) . i.e. Q m = 0.24, Q b = 0.042, S1 A = 0.76 
and h — 0.73, a normalization of as = 0.75 and a slope of the 
power spectrum of n = 0.95. The treePM-SPH code GAD GET 2 
dSDringell2005) has been applied to simulate the evolution of a cos- 
mological box with side length of Lb ox = 64fe _1 Mpc in which the 
formation a Local Group has been enforced by constraints on the 
initial conditions. There are two runs available, one with dark mat- 
ter only (DM run) and one hydrodynamical (labelled SPH run) in 
w hich we additionally fo l low th e feedback and star formation rules 
of ISnringel & Hernquistl d2003l) . as well as a uniform but e volv- 
ing ultra-violet cosmic background dHaardt &~M adau 1996]). The 
runs feature a mass resolution of tuum = 2.1 x 10 h~ M© for 
the dark matter particles (m D M = 2.54 x 1O 5 /i _1 M in the DM 
only run), m gas = 4.42 x 10 4 /i _1 Mq for the gas particles and 
wistar = 2.21 x 10 4 /i _1 Mq for the star particles. The gravita- 
tional softening length is e = O.lh^ 1 kpc = 137pc, in both the 
DM only and the SPH run. 

The stellar pop ulation synthesis model STARDUST (see 
iDevriendt et alfl999l and references therein for a detailed descrip- 
tion) has been used to derive luminosities from the stars formed 
in our simulation. This model computes the spectral energy distri- 
bution from the far-UV to the radio, for an instantaneous starburst 
of a given mass, age and metalicity. The stellar contribution to the 
total flux is calc ulated assuming a Kennicutt initial mass function 
dKennicutdl 19981) . 
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The AHF halo findeiQhas been used to identify all (sub-)haloes 
in our simulation. Note that AHF automatically (and essentially 
parameter-free) finds haloes, sub-haloes, sub-subhaloes, etc. All the 
subhaloes used in this particular study are selected in order to be 
more massive than M S ub ) 2 x 10 s h~ 1 Mg, with a peak in the 
velocity curve V ma x J> lOkm/s, and to lie within 300 kpc from 
each host's center, the latter being either the MW or M31. The 
masses of the SPH hosts, defined as the masses within a sphere 
containing A v i r ~ 390 times the cosmic mean matter density, are 
Mmw = 4.0 x 10 u /i _1 M s and M M 3i = 5.47 x 1O 11 /i~ 1 M . 
When stacking the data from the two hosts together, we found a 
total of 56 SPH and 66 DM subhaloes in this WMAP3 simulation. 
Note that our selection criterion assures that within each host a sub- 
halo contains a minimum of 1000 particles. 



related to the NFW one through p_2 = p s /4. This profile gives a 
finite total mass and its logarithmic slope decreases inwards more 
gradually than a NFW or M99 profile. When tie is large, the inner 
profile is steep and the outer profile is shallow. Typical values of 
tie found in dark matter only simulations for haloes more massive 
than 10 10 M Q are 4 < n < 7. 

Prugniel-Simien profile Finally, following the study of 
iMerritt et all d2006h . we use the an alytical approxima ti on of 
the deprojected Sersic law, given in IPrugniel & Simienl 
P&S hereafter): 



pp&s(r) 



P-2 



(3) 



3 THE DENSITY PROFILE OF SPH AND DM 
SUBHALOES 

3.1 Theoretical Models 

While it is widely accepted that a iNavarro et all l l 19961 NFW) pro- 
file provides a good description of D M haloes, it has been already 
pointed out in lDi Cintio et alj ( 1201 11) that this universal profile may 
not be the best choice when used to fit ™£>-halo densities. We will 
thus study different profiles and apply them to our simulated sub- 
structures, with particular emphasis on the density profile of sub- 
haloes in hydrodynamical simulations. 

Double-power law profiles A generalisation of the NFW profile 
is the so-called (a, ft, 7), or double power-law model: 



p a „s, 7 (V) = 



on (5-7)/ a 



(1) 



where r 3 is the scale radius and p s the scale density, characteris- 
tic of each halo a nd related to its formation time and mass (e.g . 
Prada et alj|2012l; iMunoz-Cuartas et alj|201 ll; iMaccio et ai]l2007l ; 



Bullock etalj|200ll)~ ft is a five-parameter model in which the in- 
ner and outer region have logarithmic slopes —7 and — /3, respec- 
tively, and the a parameter regulates the sharpness of the transition. 
The choice (a, /3, 7) = (1, 3, 1) provides the NFW profile, while 
(a, /5, 7) = (1.5, 3, 1.5) gives the model presented in lMoore et alj 
J1999L M99 hereafter). Besides of the NFW and M99 profiles we 
will also investigate the case of leaving the central slope as a free 
parameter, i.e. a (1, 3, 7) model. 

Einasto profile In add ition to these d ouble-power law profiles we 
test the Einasto profile I 
to the 2D Sersic model dSersicH963U ' 



( Eina stolll965 ), identical i n functional form 
; Sersidl968f ), but used instead 



to fit a spacial mass density: 



PE(r) 



P-2 



(2) 



Here r_2 is the radius where the logarithmic slope of the density 
profile equals -2 and n, also referred to as tie, is a parameter that 
describes the shape of the density profile. r_2 is equivalent to the 
scale radius r s of a NFW profile, and the density p_2 = p(r_2) is 



where n, or np&s where appropriate, is again a parameter describ- 
ing the curvature of the density profile and the quantity p is a func- 
tion of Jip&g chosen to maximize the agreement between the P&S 
model and the Sersic law. A good choice for p, when 0.6 < n < 
10, is p = 1.0 - 0.6097/n + 0.05463/n 2 jLima Neto et alj 19991) . 
used in our fitting routine. We must highlight that the shape param- 
eter riE of the Einasto profile is not the same as np&s of the P&S 
model, although they follow the same functional form. The Einasto 
profile, the P&S one and the modified NFW profile (1, 3, 7) are all 
3-parameters models. 



3.2 Application to Subhaloes 

We now apply all the above models to fit our subhaloes' density 
profiles0 The density profiles are given in radial bins logarithmi- 
cally space d from the inner rad ius compliant with the convergence 
crite rion of IPower et al. Id2003l) out to the subhaloes' edge, defined 
as m iKnollmann & Knebd fcOOsI ). The number of bins varies from 
7 for the least massive objects to 16 for the most massive ones; 
by this we assure to minimize the Poissonian noise always hav- 
ing at least 150 particles per bin, and a minimum of 1000 particles 
in total in each subhalo. We verifie d that the convergence crite- 
rion as defined in lPower et al.l ( 120031) is suitable also for subhaloes, 
and thus fully applicable to our simulation. Specifically, we used 
a lower resolution, 2048 3 particles DM-only run with three times 
higher softening length e — 411 pc, to show that the density profile 
of subhaloes converges for ~ 4.8e. This value is alway s equal or 
less than the radius found using the lPower et alj J2003l) criterion: 
all our trusted radii are thus fully converged according to the most 
conservative criterion possible and are not affected by two-body 
relaxation effects. 

We define the goodness-of-fit as 



A J 



N bin 



N bins 

2^ (logiop sim , k 
fc=l 



logiopst,k) 



(4) 



whose average value over the total number of subhaloes A 2 gives 
an indication of the fit performance. 

The results are presented in Table Q] for a WMAP3 cosmol- 
ogy, where we list the quality of fit values A 2 alongside the mean 
value of the shape parameter be, np&s or the inner (outer) slope 7 
(P) in the case of considering the double power-law models. A first 



2 http://popia.ft.uam.es/AMIGA 3 We use the IDL routine MPFIT 
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Table 1. Quality of the fits for several density profile models. The results for SPH and DM subhaloes are listed, together with the mean value of the shape 
parameter n, the inner slope 7 and the outer slope /3, respectively. 







SPH 




DM 


Profile 


A 2 


shape 


A 2 


shape 


NFW 


0.043 


7 = 1.00 


0.042 


7 = 1.00 


M99 


0.030 


7 = 1.50 


0.032 


7 = 1.50 


(1,3,7) 


0.014 


7 = 1.98 


0.015 


7 = 1.53 


(i,/U) 


0.014 


/3 = 3.80 


0.013 


/3 = 4.40 


P&S 
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Figure 1. Residuals of the density profiles of all SPH subhaloes for each of the fitted models. The mean goodness-of-fit A 2 is indicated, providing the Einasto 
model to be the best one. The radial dependence of the residuals is the same for the DM only run, thus not shown here. 



consideration regards the differences between the DM only and the 
SPH runs: we observe that the mean shape parameters in the SPH 
run are systematically higher than in the DM counterpart, which 
implies a more cuspy central slope, indicating that the net effect 
of the inclusion of baryons is a steepening of the subhaloes' den- 
sity. Our results appear to be in agreement with the pres cription of 
an adiaba tic contraction model dBlumenfhal et al]|l986h , as shown 
already in lDi Cintio et al.l d201 ll) for the most massive and most lu- 
minous subhaloes. We must remark that we are listing the average 
n over the total set o f subhaloes: there are cases, as discussed in 
iDi Cintio et al. in which the SPH subhaloes with the lowest 

baryon fraction, instead, undergo an expansion, therefore lowering 
their n. An higher shape parameter in the SPH run is also indica- 
tive of a less steep outer profile with respect to the DM only run, 
which means that tidal stripping effects are stronger on the DM 
only substructure, as reported in iLibeskind et al.l ( 2010h (see dis- 
cussion below). 

We also notice that the mean shape parameter of the P&S pro- 
file, np&s, is lower than the corresponding Einasto p arameter, riE, 
in the same run: this is expected, and found also in iMerritt et al.l 
J2005l) . A few words on the steep central slopes 7 found for the 
(1, 3, 7) model: this model imposes the outer slope to be equal to 
3, which i s not the case for subhaloes where the p rofile drops even 
faster (cf. bh et al.|[l995l; IPenarrubia et~lu1l2009l) . causing the fit- 
ting routine to provide high values of 7 when trying to adjust the 
density profile. Leaving the outer profile index f3 as a free parame- 
ter, i.e. using a (a, f3, 7) = (1, /3, 1) model, we found indeed that 
on average the outer slope of SPH subhaloes is fi = 3.8 while in 
the DM only case /3 = 4.4. This is in agreement with the average 
values obtained for the Einasto shape parameter, which is higher 



in the SPH case with respect to the DM only scenario: the higher 
riE found in the SPH subhaloes indicates a cuspy inner profile, as 
expected if adiabatic contraction is acting on the central part of the 
structures, but also an outer density profile less steep than in the 
DM case, as confirmed by the (1, /3, 1) model. 
The outer profile of subhaloes in the CLUES gas-dynamical run 
is shallower than in the pure dark matter case because of the influ- 
ence of tidal stripping whose effects, being present in both runs, are 
stronger in the DM case. Tidal stripping, which mainly acts on the 
outer part of the density profile, is able to remove more mass from 
a pure dark matter subhalo than an SPH one, owing to the deepen- 
ing of the potential in the latter case, as shown in greater detail by 
ILibeskind etall J20 lfj) . 

Finall y, we also used the expo nentially truncated profile, in- 
troduced bv lKazantzidis et al] d2004l) to deal with the divergence of 
the cumulative mass distribution of haloes as r — > 00, but we did 
not obtain improvements over the Einasto or P&S profiles. 

To further highlight the quality of the different models, in 
Fig. [T] we present the residuals between the fits and the data for 
each subhalo in the SPH simulation as a function of log(r /r_2 ) (the 
plots look akin for the DM simulation and hence are omitted). Note 
that for the (1, 3, 7) model the point where the logarithmic slope of 
the density profile equals —2 occurs at a radius r_2 = (2 — 7)7- s 
for 7 < 2: thus, for a NFW profile r_2 = r 3 and for a M99 profile 
r- 2 = r s /2. 

Neither the NFW, M99, or (1,3,7) profiles are well fitted 
over the whole radial range: while the (1, 3, 7) with a steep cen- 
tral slope may describe the data in the inner regions, it fails so in 
the outer parts. On the contrary, the radial dependence of the resid- 
uals in the Einasto model is clearly minimized with respects to all 
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Figure 2. Correlation of the shape parameter we with the subhaloes' masses in the DM only simulations, left panel, and SPH simulation, right panel. The 
error bars associated to the shape parameter are indicated for rig , as provided by the fitting routine. The lcr statistical error committed in the evaluation of rig 
is, on average, the 15% of the rig itself. The evident correlation between «e and the subhaloes' masses in the DM only as well as in the SPH run is attributed 
to the dynamical effects of tidal stripping. The additions of baryons can ulteriorly modify the density profile, specifically in its inner part. 



the other models, being consistent with zero at every radial bin. The 
case could be made that the (1, 3, 7) model performs as well as the 
Einasto model in the inner region of the density profile (which is 
also the region of interest with respect to the observations of the 
MW's dSphs). In order to assess the degree to which these results 
are affected by the choice of the radial range, we calculate the resid- 
uals, for every model, for only the innermost bins with r < r_2, 
and found that the Einasto is still the best profile overall. 

We conclude that the Einasto model outperforms all the other 
proposed profiles in terms of quality of fit, giving, over the full 
radial range, an average value of A 2 = 0.011 in the SPH run 
and A 2 = 0.012 in the DM only run. We note that the the P&S 
model also provides good results, though not as good as the Einasto 
model. 

While it is somehow obvious that these 3-parameter models 
perform better than the 2-parameter ones (such as the NFW or 
M99 model), we are reassured by the fact that even after fixing the 
shape parameter riE (np&s) of the Einasto (P&S) profile to its mean 
value, therefore reducing the free parameters to two, we still obtain 
a mean goodness-of-fit which is lower than any other 2-parameter 
model (the Einasto profile, for example, provides A 2 = 0.025 for 
the SPH run and A 2 = 0.028 for the DM run). Thus, the better 
performance of the Einasto profile is not just an artifact of having 
one free parameter more. Furthermore, our r esults are in agreeme nt 
with those of other workers in the field (e.g. lSpringel et al.l2008t) . 

In light of this we state that the need of a shape parameter n to 
fully specify the mass profile of simulated DM and SPH subhaloes 
is an indication of the non-universality of their density profiles, as 
will be highlighted in the next section. 



4 SHAPE PARAMETER - SUBHALO MASS RELATION 

In Table[T]the average values of the Einasto shape parameter tie are 
shown. However, this shape parameter varies from subhalo to sub- 
halo, spanning quite a large range 0.4 < ue< 10.4 in the SPH run 
(with a similar spread in the DM model). This naturally raises the 
question of whether this variation follows some rule or is random. 

In Fig. [2] the fitted Einasto shape parameter ns is plotted 
against subhalo mass for both the DM only run (left) and the 



SPH run (right). A clear correlation is immediately visible. This 
result forms one of the main findings of this paper: the Einasto 
shape parameter correlates directly with subhalo mass. To quan- 
tify the WE-mass correlati on, the Spearman rank coefficient S r 
dKendall & Gibbons! 1 199tj) is calculated, yielding S r = 0.70 for 
the SPH run, with a significance of practically zero confirming a 
strong correlation: the most massive objects have a higher value 
of the shape parameter, while less massive ones have smaller val- 
ues. In other words, low mass substructures are well fit by a in- 
ner density profile shallower than a NFW one, with a steep outer 
slope, while the higher mass objects are fit by a steep, cuspy-like 
inner profile. Convergence studies (wherein the number of radial 
bins used for profile fitting is drastically increased) have been per- 
formed in order to ascertain the applicability of an Einasto profile 
to our subhaloes. These tests have revealed that our fitting proce- 
dure is robust and not a result of the sampling. Similar results and 
S r values are found for the DM only run. One might argue that the 
subhalo's mass may be se en as a rather ill-d efined quantity, and a 
better proxy for mass (e.g. lKnebe et al 1201 1) should be the peak of 
the rotation curve Vm ax : when replacing the mass on the a;-axis of 
Fig. [2] with Knax we actually do not find any substantial change in 
the correlation, strongly confirming it. 

Since Fig.[2]shows the same riE-mass relation for both the DM 
only and the SPH run, we conclude that the mechanism responsible 
for this relation must be dynamical and hence is likely to be tidal 
stripping. To determine the influence of tidal stripping, the proper- 
ties of subhaloes at infall time, defined as the last time a subhalo 
crossed a sphere of physical radius 300kpc from the host's center 
0, have been examined. 

At infall time the subhaloes' density is well described by both 
an Einasto and a NFW profile, as expected for field haloes. Using 



4 The Spearman rank coefficient is a non-parametric measure of correla- 
tion that assesses how well an arbitrary monotonic function describes the 
relationship between two variables, without making any other assumptions 
about the particular nature of the relationship between the variables. The 
closer the coefficient is to 1 the stronger the correlation between the two 
variables. We use the IDL routine R CORRELATEQ to calculate it. 

5 Using the first infall time provides similar results. 
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Figure 3. Reduction of the shape parameter tie as a function of the mass 
loss between z = and Zi n f a ii for the SPH subhaloes. A similar behavior 
is found in the DM only run, thus not shown here. The best fit curve, which 
has unitary slope, is shown as a solid line. 



an Einasto profile, we find that the more mass lost since 2mfa.ii, 
the lower the value of jib at z = 0, as presented in Fig. [3] for 
the SPH case. There is thus an evident correlation between the 
amount of stripped material and the reduction of tie for each sub- 
halo. Again, this relation is quantified by the Spearman rank co- 
efficient, S r = 0.68, and showing the best curve fit, which has 
a unitary angular coefficient, as a solid line in Fig. [3] A similar 
dependence is found in the DM only run, thus not shown here, cor- 
roborating our findings that tidal stripping is the main mechanism 
able to modify the density profile of subhalo es. 



Many authors (e.g. 



l2008l : lHavashi et alj|2003l : 



Kazantzidis et alj 120041 : ISpringel etall 
Penarrubia et alfeoioh have shown that 



tidal stripping acts mainly to modify the outer region of a profile. 
Since a steepening of the outer profile entails a reduction of tie 
exactly as observed in our simulations, we therefore conclude that 
the lowering of the subhaloes' shape parameter between Zj n f a n and 
2 = is primarily due to stripping e ffects. This finding is in deed 
in agreement with the recent work of IVera-Ciro et al.l d20 1 3b who 
show that very heavily stripped objects have on average smaller tie, 
because of a steepening of the outer density p rofile. It is also worth 
noting that tidal stripping has been shown by |HavashietalJ ( l2003h 
to not only affect the subhaloes' outer regions: as a substructure 
loses mass, the central density will also decrease significantly (al- 
though the slope of the inner density profile remains unchanged). 

A word of caution is necessary at this point. The Einasto 
model's shape parameter tie describes simultaneously the slope 
of the inner and outer profile, in a single number. This can be both 
an advantage and a disadvantage. Indeed for some subhaloes in our 
simulations the value of the scale radius r_2 is close to the the in- 
nermost converged radius while the outer profile is resolved with 
many radial bins. In these cases, if the outer profile steepens yet the 
structure of the inner part remains unchanged, the Einasto fit would 
return a lower value of tie due to the fitted profile being dominated 
by the steepened outer part. This is the case, for example, of heav- 
ily stripped objects {M z=0 /M Zinfall ~ 20%) in the DM only run, 
in which the inner cuspy slope of the substructure is retained after 
infall, while the outer profile has been steepened by tidal stripping, 
thus providing a small fitted tie value. Care must be taken not to 
interpret these cases as becoming cored, in the inner region, due to 



tidal stripping. That said even in these cases, Einasto models pro- 
vide accurate fits to the density profile - more accurate than any of 
the other models as could be seen in Fig.[T]and TableQ] 

While tidal stripping is the only relevant mechanism in the 
DM only run, other effects such as baryonic feedback have to 
be taken into account in the SPH case, since they can also con- 
tribute to changes in the d ensity profile of subhaloes dZolotov et al.l 
l2012MBrooks et alj2013l) . Indeed, as opposed to the DM only sim- 
ulation, in the SPH run w e observe in some sub haloes a change 
to their inner structure (cf. |Pi Cintio et al.ll2oTH) . Since, as men- 
tioned above, tidal stripping does not change the inner slope of 
substructures, the only mechanisms able to alter the density a t 
small radii must have a baryonic origin. iDi Cintio et al. 
showed that the inclusion of baryons has a twofold effect, increas- 
ing or de creasing its central dens ity according to an adiabatic con- 
traction telumenthal et alj|l98fj|) or outflows jNavarro et"al1l 19961 : 
iGovernato et alj|2012l) model. 

To shed more light onto the effects of baryons, a one to one 
comparison of the density profiles of those subhaloes that can be 
cross-identified in th e DM and SPH run has been performed (as in 
IDi Cintio et alj20TH) . The subhaloes which experienced an expan- 
sion in the SPH run, with respect to their DM only partner, have 
lost all their gas at redshift zero, and they do not show any sign of 
star formation between Zi n f a n and z = 0. These subhaloes have an 
inner density profile shallower than the corresponding sister DM 
subhalo. On the other hand, those obje cts which have under gone 
adiabatic contraction in the SPH run (cf. lDi Cintio et alj|201lh still 
retain some gas at z — and their star formation appears to be 
on going even after infall. These subhaloes have a high value of 
the shape parameter tie, which is now well describing a steeper 
inner density profile caused by adiabatic contraction. In any case, 
the fits are still dominated by the outer profile, steepened by tidal 
stripping, where most of the bins lie and hence higher resolution 
simulations are needed to verify the effective creation of a core in 
objects with small tie values: the interplay of these two contrasting 
effects, i.e. outflows vs adiabatic contraction, and a deep analysis 
of the repercussions on the inner density profile of substructures in 
cosmological simulations will be explored in detail in a companion 
paper (Di Cintio et al., in prep.). 

Some other important conclusions can now be drawn from 
Fig. [2] Firstly, there is no evidence for any universal profile in 
simulated substructures. Secondly, we note that the majority of 
the subhaloes in both runs tend to have a small tie, while only 
the most massive ones (mostly the adiabatically contracted ones) 
have a high tie, as large as ?ie=10.4 in the SPH run and riE=9.4 
in the DM only one. This finding, as well as the goodness of the 
Einasto profile, has been con firmed from an observational point 
of view by the recent work o f lDel Popolo & Cardond I I2012I) , who 
used high quality rotation curves data of dwarf galaxies to show 
that the preferred fitting function is given exactly by the Einasto 
model and that the majority of the dwarfs tend to have shallow pro- 
files (their Fig. 3). Our mean shape parameter in the hydrodynam- 
ical run, TTe = 4.8, as well as their Einasto mean shape param- 
eter, «e,d.p. = 3.05, are both lower than the corresponding tie 
for dwarf size objects found in previou s dark matter simulations 
( iMerritt et ai1l2005l: iNavarro et al.ll2004 . In that regards we need 
to mentioned that previous results coming from such dark mat- 
ter only simulations, where cluster- and galaxy-sized haloes have 
been studi ed, showed instead a decreasing of tie for increasing 
halo m ass jNavarro et al.ll2004 IMerritt et al.ll2005t iGraham et al l 
2006bl:|Prada et alfeOOfjUGao et al.ll2008l:lHavashi & Whitell2008l: 
Navarro et al.ll2010l) . Our study indicates that there is a turnover, 
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such that the trend with mass is reversed for low mass galaxies 
(at least satellites), with both SPH and DM simulations having 
a positive correlation of he with mass. The main difference be- 
tween previous studies and this work is, besides the less massive 
objects considered here, the fact that our objects are subhaloes, and 
thus obviously affected by tidal stripping, as mentioned above. Fi- 
nally it must be noticed that the range of variation of the tie shape 
parameter found in our simulation is very large, spanning the in- 
terval 0.4 <J riE<, 10.4 in the S PH case: remarkably the sam e 
large range has also been found in iDel Popolo & Cardond d2012l) . 
with 0.29 < we. p. p . < 9.1, as well as in the recent work of 
IVera-Ciro et al.l d2013l) . based on semi-analytical mo dels of galaxy 
forma tion. Furthermore, in the observational paper of lChemin et al.l 
]]), the authors used the Einasto model to fit the rotation curves 
oftheTHiNGS^ galaxies and show that the shape parameter is near 
unity on average for intermediate and low mass halos, while it in- 
creases for higher mass haloes, being correlated with the halo virial 
mass as we find in this work. 



5 NEW OBSERVATIONAL CONSTRAINTS FOR THE 
SATELLITE GALAXIES OF THE MILKY WAY 

We now move to a practical application of our findings, only focus- 
ing on the properties of SPH subhaloes, whic h are obviously closer 
to real ity than their DM only counterparts. In lBovlan-Kolchin et al.l 
J201 l]) the observational constraints, used to establish if subhaloes 
found in cosmological simulations are possible hosts of the known 
Milky Way dwarf spheroidals, were based on the assumption that 
the underlying dark matter halo of these dSphs follows a NFW pro- 
file. Given our findings, however, it is clear that since the Einasto 
model provides the best fit to the density profile of both DM and 
SPH subhaloes, those observational constraints have to be modi- 
fiedO Note that in Eq.5-7, for clarity, we will omit the subscript e 
from the Einasto shape parameter, simply referring to it as n. The 
circular velocity of an Einasto profile follows 



v (r) oc 7(3n, a;), 



where 



7(3n,x) 



— £,3n — 1 it 

e t at 



(5) 



(6) 



is the lower incomplete gamma function and x = 2n(r /r_2) 1,/n • 
To find the radius R m ax at which dv(r)/dr — we numerically 
solve, 



-Rmax -.l/ 7l s 3n 3n- 1 / -Rmax n 3 — 2n(- 
) ) = 2 n ( ) e 

r_2 r_2 



7(3n, 2n( 

(7) 

The relation between R max and r_2, that we need in order to 
compute the observational constraints, varies depending from the 
value of the shape parameter riE (see, for example, Fig. 2 in 
lGrahametaI1l2006al) . Given the fact that t he mass density pro- 
file of the faint dSphs is sti ll uncertain (e.g. IWalker & Penarrubial 
1201 lMWolf & Bullockll2012l) we prefer to use the conservative lim- 
its given by the highest and the lowest values of the shape param- 
eter riE as obtained from our hydrodynamical simulation. In the 



6 http://www.mpi a. de/TH INGS/Overview.html 

' IVera-Ciro et al . (2013) reached similar conclusions, using semi- 
analytical galaxy formation models, while this manuscript was being pre- 
pared for submission. 



SPH run the smallest 71e= 0.4 corresponds to a relation R ma x = 
1.447r_2, while the largest ns= 10.4 gives Rmax — 2.348r— 2. 
Using these constraints, i.e. the assumption of an Einasto model 
and the corresponding range of riE-values, we computed the curves 
in the Vmax-Rmax plane for the nine brightest classical dSphs of 
the MW, namely Cvnl, Carina, Draco, Fornax, Leo I, Leo II, Sex- 
tans, Sculptor and Ursa Minor, which all have My <. —8.8 (we 
excluded Sagittarius as in iBovlan-Kolchin et alj j2o7lf) since it is 
far from dynamical equilibrium): these curves are constructed by 
normalizing each dwarf to its observationally derived values o f 
half-light mass, M1/2, and radius, ri/2, from I Wolf et al.l d2010l) . 
who showed that any uncertainty on the stellar velocity dispersion 
anisotropy is minimized at this radius, leading to accurate estima- 
tion of v(ri/2)- 

In the left panel of Fig. [4] we show the maximum circular ve- 
locity Vmax and its corresponding radius R max for all the SPH sub- 
haloes, within the MW and M3 1 hosts, whose luminosity is at least 
as high as the Draco's one, i.e. My <, —8.8. With respect to Fig. [2] 
we excluded here objects with a luminosity lower than Draco, but 
verified that the interval for the shape parameter ns is still the same. 
In Fig. [5] we also show the numerically derived Vmax-mass rela- 
tion for the SPH subhaloes, which is useful to derive the range of 
masses associated to a specific value of V m ax ■ The grey symbols 
in Fig. [4] correspond to the subhaloes that are brighter than For- 
nax, which is the brightest classical dwarf considered here to con- 
struct the observational constraints having My = —13.2, and the 
black circles indicate all the remaining subhaloes with luminosity 
— 13.2 < My < —8.8. We also plot the newly constrained obser- 
vational limits, as solid lines, coming from the assumption that the 
MWs dSphs are embedded in haloes that follow the Einasto profile 
with varying shape parameter tie between 0.4 < tie < 10.4, and, 
as dashed lines, the previously used constraints coming from the 
NFW model. 

We observe that, while the employment of an Einasto profile 
leads to a good agreement between observations and the SPH sub- 
haloes0 it still appears to be not sufficient to explain the R m ax- 
V^nax pairs of the most massive SPH subhaloes which still lie in 
the lower right part of the plane, outside the constraints. Note that 
with a shape parameter varying between 0.4 < he < 10.4 we have 
allowed the observational constraints to cover a wider range, in the 
Rmax- Vmax plane, with respect to the NFW constraints, but even 
this assumption is not enough to reconcile simulation and observa- 
tion. However, those massive SPH objects, which we color-coded 
in grey, appear to have a luminosity, My < —13.2, not compati- 
ble with any of the satellites used to derive the observational con- 
straints. 

We remind the reader that the data plotted in Fig. prefer to 
both the MW and M31 galaxies: a total number of 14 subhaloes 
brighter than Fornax is thus found within the two hosts (and only 1 1 
if we relax the magnitude cut from My < —13.2 to My < —14). 
For each host halo we therefore have 5 to 7 objects brighter than 
the classical dwarfs used to compute the observational constraints. 
Three of them may be associated with the Large and Small Mag- 
ellanic Clouds (LMC and SMC, respectively) and the Sagittarius 
galaxy. In fact, objects with My <^ —16.2 can be conservatively 
considered as the analogous of the LMC and SMC: we should 
therefore exclude these simulated subhaloes from the discussion. 

With these associations there are, within each host, only two 



8 We find a complete agreement between observations and DM only sub- 
haloes, which are though not shown in the plot for clarity. 
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Figure 4. Left panel, Vmax-fimax pairs for the SPH subhaloes within the MW and M31 hosts. The subhaloes have been color-coded by their luminosity: in 
black the ones corresponding to the luminosity of the observed classical dSphs, in grey the ones which are brighter than My = — 13.2. The 2cr observational 
constraints for the MW's dSphs are indicated as a solid line for the Einasto profile with shape parameter 0.4 < tie < 10.4, obtained considering all the 
subhaloes with My < —8.8, and as dashed line for the NFW profile. Right panel, same criterion to color-code the subhaloes, although this time the 2<r 
observational constraints for the MW's dSphs are derived using an Einasto profile with shape parameter 1.6 < tie < 5.3, obtained considering only the 
subhaloes with -13.2 < My < -8.8. 



to four subhaloes left with —16.2 <C My <C —13.2 which do not 
have a counterpart in the real universe: such a small sample, i.e. 
the ~ 10% over the total number of objects found within each halo 
in the SPH run, can be explained as a statistical fluctuation due 
to our small number statistics. In the future, to confirm this halo- 
to-halo variation in the subhaloes population, it will be necessary 
to study many realization of a high resolution MW-like object. In 
this work the luminosity function, averaged over the MW and M31 
subhaloes, has been shown to be in a greement with the observa - 
tional data of Milky Way like galaxies dStrigari & Wechslei1l2012l) . 
while slightly deviati ng from the Milky Way itself in the interval 
-16 < My < -13 dKnebe et alj|201lh exactly because of these 
two to four overabundant objects in this range. Moreover, we re- 
mind that our simulation also reproduces the luminosity vs veloc- 
ity disper sion correlation obs erved for the satellite galaxies of MW 
and M31 dWalker et al]|2009h , as shown in lKnebe et aT] feoill) . 

As an additional remark, we note that most of the brightest 
subhaloes in our simulation are the ones that experienced adiabatic 
contraction, being situated in the lower right part of the Vmax-fimax 
lane at redshift z = 0, as studied and explained in lDi Cintio et al .1 
201 lh . These subhal oes are substant i ally d ifferent from the ones 
found in the work of IVera-Ciro et al.l d201 3b . who used dark mat- 
ter only simulations with semi-analytical galaxy formation models 
that do not show adiabatic contraction: their most luminous, bright- 
est objects are found in the upper-right of the V ma x-flmax plane, 
contrary to what we obtained in our hydrodynamical simulations. 

We now proceed to again compute the shape parameter range 
based only upon those subhaloes that satisfy the luminosity re- 
quirement, i.e. those objects whose My is within the range of the 
observed dSphs luminosity. We were therefore able to restrict the 
range, finding a shape parameter lying within 1.6 < he < 5.3, 
with a mean value We = 3.2. In the right panel of Fig.|4]we use the 
same black-grey colouring scheme for the SPH subhaloes as be- 
fore, and we plot the observational limits based on the newly con- 
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Figure 5. Correlation of the subhaloes maximum circular velocity with the 
subhaloes' masses in the SPH simulation. 



strained range for the shape parameter 1.6 < tie < 5.3. The result 
is a perfect agreement between the expected Kn ax -7? m ax values of 
the observed dSphs and the V max --Rmax pairs of the simulated sub- 
haloes with corresponding luminosities. 

We conclude that our findings, based upon self -consistent hy- 
drodynamical simulations of a constrained Local Group in a cos- 
mological context, strongly supports the notion that the observed 
satellite galaxies of the Milky Way are actually compatible with 
being embedded in dark matter haloes whose density profiles show 
considerable differences, following an Einasto model with 1.6 < 
riE < 5.3 and mean value tie = 3.2: the majority of the dSphs 
may have an inner profile shallower than the previously assumed 
NFW one while an outer profile steepened by tides. 
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6 CONCLUSION AND DISCUSSION 

Using a constrained simulation of the Local Group, performed 
within the CLUES project, it has been shown that: 

• the density profile of subhaloes in both dark matter only and 
hydrodynamical simulations is best approximated by an Einasto 
profile in which the shape parameter tie is free to vary, and that 

• there is a clear trend of growing tie with increasing subhaloes 
mass, in both the dark matter only and hydrodynamical run. 

The structural effect associated with tidal stripping is likely the 
main mechanism able to modify the subhaloes' density profile: the 
effect of mass loss due to tidal stripping is the reduction of the shape 
parameter tie between the infall and the present time. A correlation 
betwee n he and the amount o f stripped material has also been ar- 
gued in lVera-Ciro et al. Differences in the inner profile of 
subhaloes, between the pure DM and SPH run, can instead be at- 
tributed to baryonic processes dDi Cintio et al Hoi ]]): they result in 
adiabatic contraction of the dark matter halo when gas is retained 
in the central regions and star formation is still on going after in- 
fall, and in expansion when gas is removed from the dwarfs due to 
stellar feedback and ram pressure stripping, with no signs of star 
formation after the subhalo has entered the host's halo. These bary- 
onic effects, acting mainly on the inner part of the density profile, 
do not modify the overall tie -mass relation which is driven by tidal 
stripping also in the SPH case. 

The majority of our SPH subhaloes have a small 
we, as reported in the rig ht panel of Fig. [2] remarkably, 
Del Popolo & Cardond J2012h found similarly small values of tie 
in observed dwarf galaxies, using high quality rotation curves. 
Moreover, evidences of the fact that at least some of the MW's 
dSphs may have a shallow profile, co mpatible with a small tie, are 
given in lWalker & Penarrubial d201 ll) . who showed that the profiles 
of the Fornax and Sculptor dSphs are consistent with cores of 
constant density at a high confidence level. Nevert heless, the actual 
mass p rofile of the MW's dSphs is still uncertain: IWolf & Bullockl 
horn claimed that, even with an isotropic velocity dispersion, not 
all the dSphs prefer constant-density cores and that, instead, some 
of them favor a cuspy inner profile. All these findings do no longer 
support the notion of a universal subhalo mass profile; subhaloes 
of differing mass cannot be rescaled to have self-similar profiles: 

their mass (or so to speak size) matters. 

In light of these results we revisited the lBovlan-Kolchin et al .1 

d201ll) observational limits for possible hosts of the MW's dSphs, 
assuming that the latter are embedded in haloes that follow an 
Einasto profile, as opposed to the earlier assumption of NFW pro- 
file, with variable shape parameter tie, and using the conserva- 
tive limits 0.4 <J tie <C 10.4 provided by our hydrodynamical 
simulations. While using the Einasto profile is enough to com- 
pletely explain the maximum velocity of the most massive DM 
subhaloes, an issue still remains with respect to the most mas- 
sive SPH subhaloes: these objects experienced adiabatic contrac- 
tion dDi Cintio et al.|[2oTTh and their -Rmax-Vmax pairs are still ly- 
ing outside the expected observational constraints. However, these 
subhaloes appear to be brighter than Fornax, which is the brightest 
dSph used when constructing the observational constraints: thus, 
they should not be considered in the comparison. Once the Large 
Magellanic Cloud, Small Magellanic Cloud and Sagittarius galaxy 
analogues are removed, we still have two to four unaccounted ob- 
jects per halo, whose luminosity is higher than the luminosity of 
the classical dwarfs: we argue that, being only 10% of the total set, 
they can be interpreted as a statistical deviation. 



Leaving only the SPH subhaloes with -13.2 < M v < -8.8, 
i.e. those in agreement with the luminosity of the nine classi- 
cal dSphs, we show that an Einasto profile with shape parameter 
1.6 < tie < 5.3 provides an accurate matching between simula- 
tions and observ ations, alleviating the "massi ve failures" problem 
first addressed in lBovlan-Kolchin et alj d20 1 lh - The mean value of 
the shape parameter for them is tie = 3.2, indicating that the ma- 
jority of the MW's satellite galaxies are consistent with dark matter 
haloes whose profile is an Einasto one, steepened outside by the 
effects of tidal stripping and possibly shallower than the previously 
accepted NFW towards the center. 



We further note that our simulated host haloes masses are 
at the low end of current observational estimates, and this may 
be one of the reason for having only a few objects i n the lumi- 
nosity range -16 < My < -13: |Pi Cintio et al.1 ( feoill) sug- 
gested that the host halo mass is directly connected to the num- 
ber of massive subhaloes found in simulations, when comparing 
the results of dark matter only simulations based on a WMAP3 
versus WMAP5 cosmology, the latter showing a higher host halo 
mass and consequently a higher number of massive subhaloes. 
The dependence of the number of "too massive su bhaloes" on 
halo mass has been further explored and quantified by IWang et all 
d2012h . lBovl an-Kolchin et alj (201 ll ) extensively discussed the pos- 
sibility that the reduction of the Milky Way m ass could solve 
the problem, and recently IVera-Ciro et al.l d2013l) suggested that 
a Milky Way mass ~ 8 ■ 10 11 Mq provides a good match be- 
tween observations and semi-analytical galaxy formation mod- 
els. In our simulations we have slightly lower masses for the 
MW and M31, between 5.5 and 7.5 ■ 10 n Me, these values be- 
ing at th e low end of mass estimates obt a ined using different 



20101: 



methods l Karachentsev & Kashib adze 200fj; IWatkins et al 
IPeason et al J2012h . According to the model of lWang et al.l 1 2012b 
a MW mass of M M w = 5.5 ■ IO^Mq will give a 84% probabil- 
ity of finding only three satellite galaxies with a circular velocity 
peak higher than Umax > 30 km/s: this is what we expect since, 
apart from the LMC, SMC a nd Sagittarius, all t he other classical 
dwarfs have been shown (e.g. lStrigari et al.l2010h to inhabit haloes 
with a maximum circular velocity below 30 km/s. However, such 
a low mass for the Milky Way will reduce the probability that it 
hosts two satellites as the LMC and SMC; that our galaxy system 
is rare, with only ~ 3.5% of the MW-like candidates having two 
satellites as bright as the Magellanic Clouds, has been fou nd obser- 



vationally by studies using the Sloan Pigital Sky Survey ( Liu et al.1 
|2011| ; iGuo et al.ll20lll ; lLares et al.ll201ll : lTollerud et al.ll201 ll) . We 
remark that, despite a lower Milky Way mass, subhalo density 
profiles should nevertheless be described by an Einasto model in 
order to properly match the kinematic of the observed classical 
dSphs with the subhaloes in hydrodynamical simulations. Finally, 
even assuming a small mass for the Milky Way and an Einasto 
profile for its satellite galaxies, there is still a problem in assign- 
ing the correct halo masses to dwarf galaxies, as highlighted in 
iFerrero et al.l d2012l) . These authors showed that the MW's dSphs, 
as well as many isolated dwarf galaxies with spatially resolved rota- 
tion curves and stellar mass 10 6 < M ga i/MQ < 10 7 , seem to live 
in haloes with M < 10 10 M(T), which is a t odd s with the abundance - 
matching prediction of IGuo et all d2010h and lMoster et al.l d2oToh : 
the validity of ACDM at such scales is still disputable. 
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